% "Repricing Avalanches" Online Appendix D, Model with constant money
% growth policy.
% Baseline code "simulations.m" is modified to allow for variable real
% wage. See Appendix D for explanation.
% Modified parts are marked by "%----".
% Need subroutine get_equilibrium.m and setparams.m
% Choose gamma_m in line 22 and a saved file name at the last line.
% October 2022
% Makoto Nirei

% State variable is modified to r = p/w
% Read pb as w * \undreline{r} and pt as w * r^*.
% Then, the computation of pb and pt is unalterned.
% Keep track of w*r when w varies as m varies.

clearvars;

% Set parameter values by setparams.m
setparams;

% three calibrations for money demand elasticity
%gamma_m=0.0;
gamma_m=0.1;
%gamma_m=0.5;

piv=0.03; %[0.01 0.03 0.05 0.07 0.09 0.12];   % inflation levels
pilength=length(piv);

% Calvo timings. The same random series is reused for various levels of pi
Fcalvo=rand(Z,1);
calvot = log( 1-Fcalvo )/(-mu*n); % Elapsed time between two Calvo shocks
calvot = cumsum(calvot); % Calendar time of Calvo shocks
maxmonth=floor(calvot(Z)/month_time);   % # of months in the simulation

% Calvo and productivity shocks. Reused for various pi's
ind_c(:,1) = ceil(rand(Z,1)*n); % index of Calvo-hit firms
ind_c(:,2) = (rand(Z,1)<zeta); % productivity shock. 0=non-switch, 1=switch

calvom=zeros(maxmonth,1);   % calvom(i): cumulative # of Calvo shocks by i-th month
calvom(1)=sum(calvot<month_time); % calvom(1)= # Calvo in the 1st month
for i=2:maxmonth
    calvom(i) = sum(calvot<(month_time*i));
end

% main loop for each inflation level piv(i_pi)
dlogP=zeros(Z,pilength); nump=zeros(Z,pilength); NpY=zeros(Z,pilength); % preallocation
negnump=zeros(1,pilength); ss_pi=cell(pilength,1);
for i_pi=1:pilength
    tic;
    pi=piv(i_pi);

    get_equilibrium; % Compute steady state
    ss_pi{i_pi}=ss;

    disp('eta rho zeta mu delta w');
    disp([eta rho zeta mu delta w]);
    disp('a(i) pbar(i) pstar(i) q(i) pu(i)'); % a, \underline{p}_a, p^*_a, \bar{p}_a
    disp([a(1) pb(1) pt(1) q(1) pu(1); a(2) pb(2) pt(2) q(2) pu(2)]);
    disp(pi);

    % stationary distribution g(s), using varsigma=0.5
    Gs = @(s,mu,pi,qL,qH) ((qH.^(s*mu/pi)-1)/(qH^(mu/pi)-1) + (qL.^(s*mu/pi)-1)/(qL^(mu/pi)-1))/2;

    % shorthands
    ometa=1-eta; % one-minus-eta
    obometa=1/ometa;
    epb=(pb').^ometa; % "ep" = p^{1-\eta}. "pb"=p underbar. "pt"=p top=p*.
    ept= (pt').^ometa;
    eq=(q').^ometa;
    epu=(pu').^ometa; % "pu"=\bar{p}
    equ=(qu').^ometa; % "eq"=q^{1-\eta}
    l_pb=log(pb);

    psimu_tmp=zeros(n/2,2); psimu=zeros(n,2); epsimu=zeros(n,2);
    % draw initial price profile for each a(i)
    % Using G(s|a) = (q_a^{(mu/pi)s}-1)/((mu/pi)log(q_a))
    % Letting u=rand=G(s|a), obtain p= pb_a q_a^s
    for i=1:2
        psimu_tmp(:,i)=((rand(n/2,1))*(q(i)^(mu/pi)-1)+1).^(pi/mu)*pb(i);
    end
    psimu_tmp = psimu_tmp / mean(mean(psimu_tmp.^(1-eta)))^(1/(1-eta)); % normalize to P=1
    psimu_init=psimu_tmp(:);
    psimu(:,1)=psimu_init;
    epsimu(:,1)=psimu(:,1).^ometa;  % epsimu(:,1) = p_i^{1-\eta}
    epsimu(1:n/2,2)=1;              % epsimu(:,2) = a_i
    epsimu(n/2+1:end,2)=2;          % Initially, the first half firms are a_L, the second half a_H.
    epsold=epsimu(:,1);

    %----NEW: introduce new variables
    prev_jump_time = 0; % time of last price jump
    [gap_to_floor, ind_gap_to_floor]=min(log(epb(epsimu(:,2))) - log(epsimu(:,1)));
    if gap_to_floor<0, disp('error: gap_to_floor<0'); pause; end
    dlogP_sub = [];
    %----

    for z=1:Z % z-th Calvo shock
        disp(z);

        %----NEW: adjustment by real wage
        shift_r_by_w = (calvot(z)-prev_jump_time)*pi*gamma_m*(eta-1); % increase of log(r^{1-\eta}) by increase of w if no jump occurs before next Calvo
        while  shift_r_by_w > gap_to_floor && shift_r_by_w>0
            % price floor hit by wage decline
            prev_jump_time=prev_jump_time + gap_to_floor/pi/gamma_m/(eta-1);
            epsimu(:,1) = epsimu(:,1) * exp( gap_to_floor );

            ePold = mean(epsimu(:,1));  % aggregate P^{1-eta}(t-)
            epsold(:,1) = epsimu(:,1);
            cfirm=ind_gap_to_floor;
            epsimu(cfirm,1) = ept(epsimu(cfirm,2)); % reprice to p^*(a)
            ePnew = mean(epsimu(:,1));  % update eP = P^{1-eta}
            deP = ePnew/ePold;      % shift in eP = P^{1-\eta} from eP(-)
            ind2 = find(epsimu(:,1)/deP > epb(epsimu(:,2))); % update ep and find firms hitting threshold.
            % note that eP shifts down when eP shifts up since 1-eta<0
            numf = 1;   % number of repricing firms including the triggering firm
            epsimu(ind2,1)=epsimu(ind2,1).*eq(epsimu(ind2,2)); % "ind2" firms reprice to pq
            flag = length(ind2);    % # of ind2 firms
            while flag>0 && numf < n    % Avalanche continues until no ind2 firms
                numf=numf+flag;         % accumulate # of repricing firms
                ePnew = mean(epsimu(:,1));  % Update eP
                deP_tmp=ePnew/ePold;        % shift in eP from eP(t-)
                ind2 = find(epsimu(:,1) > deP_tmp * epb(epsimu(:,2)));  % find the next round ind2
                epsimu(ind2,1) = epsimu(ind2,1).*eq(epsimu(ind2,2));    % update ep
                flag=length(ind2);      % # of ind2 firms
            end
            deP = ePnew/ePold;
            epsimu(:,1) = epsimu(:,1)/deP;      % Update the price profile

            %----NEW: adjust by real wage: r^i = (P^i/P)/w. w declines as m since P rises
            epsimu(:,1) = epsimu(:,1) * deP^gamma_m;
            %----

            dlogP_sub = [dlogP_sub;
                [z log(deP)*obometa numf shift_r_by_w gap_to_floor prev_jump_time]];   % Store the price increase at this instant
            [gap_to_floor, ind_gap_to_floor]=min(log(epb(epsimu(:,2))) - log(epsimu(:,1)));
            shift_r_by_w = (calvot(z)-prev_jump_time)*pi*gamma_m * (eta-1);
        end
        epsimu(:,1) = epsimu(:,1) * exp( (calvot(z)-prev_jump_time)*pi*gamma_m*(eta-1) );
        %----

        ePold = mean(epsimu(:,1));  % aggregate P^{1-eta}(t-)
        epsold(:,1) = epsimu(:,1);

        cfirm=ind_c(z,1);           % index of z-th Calvo-hit firm
        if ind_c(z,2)==0            % No switch of productivity
            epsimu(cfirm,1) = ept(epsimu(cfirm,2)); % reprice to p^*(a)
            ePnew = mean(epsimu(:,1));  % update eP = P^{1-eta}
            deP = ePnew/ePold;      % shift in eP = P^{1-\eta} from eP(-)
            ind2 = find(epsimu(:,1)/deP > epb(epsimu(:,2))); % update ep and find firms hitting threshold.
            % note that eP shifts down when eP shifts up since 1-eta<0
            numf = 1;   % number of repricing firms including Calvo
            epsimu(ind2,1)=epsimu(ind2,1).*eq(epsimu(ind2,2)); % "ind2" firms reprice to pq
            flag = length(ind2);    % # of ind2 firms
            while flag>0 && numf < n    % Avalanche continues until no ind2 firms
                numf=numf+flag;         % accumulate # of repricing firms
                ePnew = mean(epsimu(:,1));  % Update eP
                deP_tmp=ePnew/ePold;        % shift in eP from eP(t-)
                ind2 = find(epsimu(:,1) > deP_tmp * epb(epsimu(:,2)));  % find the next round ind2
                epsimu(ind2,1) = epsimu(ind2,1).*eq(epsimu(ind2,2));    % update ep
                flag=length(ind2);      % # of ind2 firms
            end
        else    % productivity switch
            if epsimu(cfirm,2)==1    % Shock hits a_L firm. Positive tech shock; price cut
                epsimu(cfirm,2)=2;  % Update to a_H
                epsimu(cfirm,1)=ept(2); % Reprice to p^*(a_H)
                ePnew = mean(epsimu(:,1));
                numf = -1;          % Negative avalanche = price cuts
                deP = ePnew/ePold;
                ind2 = find(epsimu(:,1)/deP < epu(epsimu(:,2)));    % Upper threshold for p; Lower threshold for ep
                epsimu(ind2,1) = epsimu(ind2,1).*equ(epsimu(ind2,2));   % Update p using pu
                flag = length(ind2);                                % number of price cutting firms
                while flag ~= 0 && -numf < n
                    numf=numf - flag;                               % price cut is stored as a negative numf
                    ePnew = mean(epsimu(:,1));
                    deP_tmp=ePnew/ePold;
                    ind2 = find(epsimu(:,1) < deP_tmp * epu(epsimu(:,2)));  % using threshold pu
                    flag=length(ind2);
                    epsimu(ind2,1) = epsimu(ind2,1).*equ(epsimu(ind2,2));
                end
            else                    % Shock hits a_H firm. Negative shock. Price increase.
                epsimu(cfirm,2)=1;  % Update to a_L
                epsimu(cfirm,1)=ept(1); % Reprice to p^*(a_L)
                ePnew = mean(epsimu(:,1));
                numf = 1;
                deP = ePnew/ePold;
                ind2= find(epsimu(:,1)/deP > epb(epsimu(:,2)));
                epsimu(ind2,1)=epsimu(ind2,1).*eq(epsimu(ind2,2));
                flag=length(ind2);
                while flag>0 && numf < n
                    numf=numf+flag;
                    ePnew = mean(epsimu(:,1));
                    deP_tmp= ePnew/ePold;
                    ind2 = find(epsimu(:,1) > deP_tmp * epb(epsimu(:,2)));
                    epsimu(ind2,1) = epsimu(ind2,1).*eq(epsimu(ind2,2));
                    flag=length(ind2);
                end
            end
        end
        deP = ePnew/ePold;
        epsimu(:,1) = epsimu(:,1)/deP;      % Update the price profile for the next Calvo event

        %----NEW: adjust by real wage
        epsimu(:,1) = epsimu(:,1) * deP^gamma_m;
        %----

        dlogP(z,i_pi) = log(deP)*obometa;   % Store the price increase at this instant
        nump(z,i_pi)=numf;                  % Store # of repciring firms at this instant including Calvo-hit firm
        NpY(z,i_pi)= mean(epsimu(:,1).^(-eta/(1-eta)) ./a(epsimu(:,2))' );  % N/Y

        %----NEW: adjust by real wage, and update time
        [gap_to_floor, ind_gap_to_floor]=min(log(epb(epsimu(:,2))) - log(epsimu(:,1)));
        prev_jump_time = calvot(z);
        %----
    end
    negnump(i_pi)=sum(nump(nump(:,i_pi)<0,i_pi))/sum(abs(nump(:,i_pi)));    % Fraction of price cuts
    toc;
end
%save('variable_wage_01.mat','-v7.3');